Purpose

This document will a modified workflow from my work on the India Emissions Project.

Includes creating buffers around Urban areas at specified steps. Start with 1km as proof of concept on the subset data.

Load required packages. Mapview and Mapedit are data visualization packages for exploring the data. There are some simple maps at the end of the script that use these packages. They are not essential for analysis.

require(sp)
require(sf)
require(data.table)
require(tidyverse)
require(mapview)
#require(mapedit)

Read in the file of india and filter into urban and rural areas. Data originally came in R native format. Subset used for testing purposes After initially exploring the data in ArcMap I exported a subset as a shapefile. A similar workflow could be accomplished by subsetting based on a geographic state of India.

# Subset created from the shapefile and used for testing
sub<-st_read("india/subset", layer = "india_subset")
## Reading layer `india_subset' from data source `/Users/maxblasdel/Desktop/India Emissions/india/subset' using driver `ESRI Shapefile'
## Simple feature collection with 1001 features and 58 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 77.84991 ymin: 19.06007 xmax: 79.97593 ymax: 19.91705
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
urb<-sub %>%
  filter(TRU == "Urban")
## Warning: package 'bindrcpp' was built under R version 3.4.4
rur<-sub %>%
  filter(TRU == "Rural")

Initial explore of the data

sp::plot(sub['TRU'], main = 'Urban Rual Distinction in India Subset')

The data needs to be reprojectd so that meters buffer can be applied.

Projected to Asia South Albers Equal Area Conic EPSG: 102028

urb<-st_transform(urb, crs = 102028)
rur<-st_transform(rur, crs = 102028)

Create buffers of different radius. Units are in meters. The union function combines overlapping buffers from urban areas that are close together. This is important for the western urban area which is actually two urban polygons that border each other. Without the st_union() function the st_buffer() would create some odd overlaps that cause problems for analysis.

# Define buffer distance.
buffer_distances<-as.list(c(1000, 2000, 3000, 4000, 5000))

bufs <- lapply(buffer_distances, function(x) {
  st_buffer(urb$geometry, x) %>% 
    st_union()
})

Add a new field to the rural polygons showing their area. This will be used for analysis with the buffers by calculating how much overlap there is from the urban buffers.

rur<-rur %>%
  mutate(area = st_area(rur)) # %>% as.numeric())

Look into how these intersect. I need to determine the area of the intersecting polygons from the buffer. Attribute values are kept constant and need to be scaled based on the new area of the polygons.

intersections <- lapply(bufs, function(x) {
  st_intersection(rur, x)
})

#Plot check
mapview(intersections[[5]])

Calculate the area of the buffered polygon intersection.

intersections_area <- lapply(intersections, function(x) {
  x %>% 
  mutate(new_area = st_area(.)) %>% #as.numeric()) %>% 
  mutate(proportional_area = new_area/area)
})

# Check outputs
intersections_area[[5]][,59:61] %>% head()
## Simple feature collection with 6 features and 3 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -5337753 ymin: 3437207 xmax: -5327275 ymax: 3441279
## epsg (SRID):    102028
## proj4string:    +proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
##            area         new_area proportional_area
## 1 8536466 [m^2]   12371.47 [m^2]    0.00144925 [1]
## 2 4181315 [m^2] 1911339.43 [m^2]    0.45711440 [1]
## 3 5867574 [m^2] 4176486.77 [m^2]    0.71179104 [1]
## 4 4248200 [m^2] 1983371.32 [m^2]    0.46687337 [1]
## 5 6473603 [m^2] 6473594.10 [m^2]    0.99999856 [1]
## 6 3596672 [m^2]  627233.72 [m^2]    0.17439280 [1]
##                         geometry
## 1 POLYGON ((-5334713 3440767,...
## 2 POLYGON ((-5330057 3441023,...
## 3 POLYGON ((-5332211 3441263,...
## 4 POLYGON ((-5335173 3440697,...
## 5 POLYGON ((-5335164 3440698,...
## 6 POLYGON ((-5327275 3438977,...

Potential workflows with new calculations.

We have just calculated buffer areas around polygons that have been identified by a given attribute. These buffers are essentially geographic areas of influence centered on urban areas. We can use these buffers to calculate proportions of population that will be potentially affected by an action. In this project we are concerned with the number of people in the rual polygons that are within each buffer. This method assumes equally distributed population through each rural polygon, which is often the best you can do with census level polygons.

Each variable of interest can be multiplied by the proportional area to produce a new proportionally weighted variable. In this example eg_ck_r is a measure of per capita energy use.

intersections_area[[1]] %>% 
  transmute(eg_ck_r_new = eg_ck_r*proportional_area,
            eg_ck_r = eg_ck_r,
            area = new_area,
            proportional_area = proportional_area) %>% 
  head()
## Simple feature collection with 6 features and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -5336528 ymin: 3433379 xmax: -5329169 ymax: 3437279
## epsg (SRID):    102028
## proj4string:    +proj=aea +lat_1=7 +lat_2=-32 +lat_0=-15 +lon_0=125 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
##     eg_ck_r_new eg_ck_r            area proportional_area
## 1 0.3476676 [1]   12.95  209964.6 [m^2]    0.02684692 [1]
## 2 7.0348769 [1]   12.95 4047812.4 [m^2]    0.54323375 [1]
## 3 6.0909571 [1]   12.95  985452.8 [m^2]    0.47034419 [1]
## 4 2.3495462 [1]   12.95  784806.9 [m^2]    0.18143214 [1]
## 5 5.7139737 [1]   12.95 1428957.5 [m^2]    0.44123350 [1]
## 6 3.0587506 [1]   12.95 3217312.8 [m^2]    0.23619696 [1]
##                         geometry
## 1 POLYGON ((-5331938 3437274,...
## 2 POLYGON ((-5332571 3436961,...
## 3 POLYGON ((-5330764 3436951,...
## 4 POLYGON ((-5330287 3435854,...
## 5 POLYGON ((-5335528 3436138,...
## 6 POLYGON ((-5331264 3435586,...

Writing the intersection outputs out for later use. These can be connected to the polygons by the C_CODE0 column. This might want to be done since the operations in this markdown can take a while to complete.

# st_write(intersections_area[[1]], "outputs", layer = "intersection_one.shp", driver = "ESRI Shapefile")

Mapping some of the results. Mapping Check a fairly heavy lift for the full dataset (crashes laptop)

# Define colors
cols<-sf.colors(n=5, categorical = T)

mapview(bufs[[5]], col.regions=cols[1]) +
mapview(bufs[[4]], col.regions=cols[2]) +
mapview(bufs[[3]], col.regions=cols[3]) +
mapview(bufs[[2]], col.regions=cols[4]) +
mapview(bufs[[1]], col.regions=cols[5]) +
  mapview(rur, legend=F)

Mapview package has some nice features for displaying data. Here we are syncing two maps together to better see the area included in the buffers compared with the rural polygon file.

m1<-mapView(intersections[[5]], legend=F, burst=F)
m2<-mapview(rur, legend = F)

sync(m1,m2)